home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / siman / siman.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  6.9 KB  |  243 lines
  1. /* siman/siman.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdio.h>
  22. #include <math.h>
  23. #include <stdlib.h>
  24. #include <string.h>
  25. #include <assert.h>
  26.  
  27. #include <gsl/gsl_rng.h>
  28. #include <gsl/gsl_siman.h>
  29.  
  30. /* implementation of a basic simulated annealing algorithm */
  31.  
  32. void 
  33. gsl_siman_solve (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
  34.          gsl_siman_step_t take_step,
  35.          gsl_siman_metric_t distance,
  36.          gsl_siman_print_t print_position,
  37.          gsl_siman_copy_t copyfunc,
  38.          gsl_siman_copy_construct_t copy_constructor,
  39.          gsl_siman_destroy_t destructor,
  40.          size_t element_size,
  41.          gsl_siman_params_t params)
  42. {
  43.   void *x, *new_x;
  44.   double E, new_E;
  45.   int i, done;
  46.   double T;
  47.   int n_evals = 0, n_iter = 0, n_accepts, n_rejects, n_eless;
  48.  
  49.   /* this function requires that either the dynamic functions (copy,
  50.      copy_constructor and destrcutor) are passed, or that an element
  51.      size is given */
  52.   assert((copyfunc != NULL && copy_constructor != NULL && destructor != NULL)
  53.      || (element_size != 0));
  54.  
  55.   distance = 0 ; /* This parameter is not currently used */
  56.  
  57.   if (copyfunc) {
  58.     x = copy_constructor(x0_p);
  59.     new_x = copy_constructor(x0_p);
  60.   } else {
  61.     x = (void *) malloc (element_size);
  62.     memcpy (x, x0_p, element_size);
  63.     new_x = (void *) malloc (element_size);
  64.   }
  65.  
  66.   T = params.t_initial;
  67.   done = 0;
  68.  
  69.   if (print_position) {
  70.     printf ("#-iter  #-evals   temperature     position   energy\n");
  71.   }
  72.  
  73.   while (!done) {
  74.     E = Ef (x);
  75.     n_accepts = 0;
  76.     n_rejects = 0;
  77.     n_eless = 0;
  78.     for (i = 0; i < params.n_tries - 1; ++i) {
  79.       if (copyfunc) {
  80.     copyfunc(x, new_x);
  81.       } else {
  82.     memcpy (new_x, x, element_size);
  83.       }
  84.  
  85.       take_step (r, new_x, params.step_size);
  86.       new_E = Ef (new_x);
  87.       ++n_evals;        /* keep track of Ef() evaluations */
  88.       /* now take the crucial step: see if the new point is accepted
  89.      or not, as determined by the boltzman probability */
  90.       if (new_E < E) {
  91.     /* yay! take a step */
  92.     if (copyfunc) {
  93.       copyfunc(new_x, x);
  94.     } else {
  95.       memcpy (x, new_x, element_size);
  96.     }
  97.     E = new_E;
  98.     ++n_eless;
  99.       } else if (exp (-(E - new_E)/(params.k * T))
  100.          * gsl_rng_uniform(r) < 0.5) {
  101.     /* yay! take a step */
  102.     if (copyfunc) {
  103.       copyfunc(new_x, x);
  104.     } else {
  105.       memcpy(x, new_x, element_size);
  106.     }
  107.     E = new_E;
  108.     ++n_accepts;
  109.       } else {
  110.     ++n_rejects;
  111.       }
  112.     }
  113.  
  114.     if (print_position) {
  115.       /* see if we need to print stuff as we go */
  116.       /*       printf("%5d %12g %5d %3d %3d %3d", n_iter, T, n_evals, */
  117.       /*           100*n_eless/n_steps, 100*n_accepts/n_steps, */
  118.       /*           100*n_rejects/n_steps); */
  119.       printf ("%5d   %7d  %12g", n_iter, n_evals, T);
  120.       print_position (x);
  121.       printf ("  %12g\n", E);
  122.     }
  123.  
  124.     /* apply the cooling schedule to the temperature */
  125.     /* FIXME: I should also introduce a cooling schedule for the n_tries */
  126.     T /= params.mu_t;
  127.     ++n_iter;
  128.     if (T < params.t_min) {
  129.       done = 1;
  130.     }
  131.   }
  132.  
  133.   /* at the end, copy the result onto the initial point, so we pass it
  134.      back to the caller */
  135.   if (copyfunc) {
  136.     copyfunc(x, x0_p);
  137.   } else {
  138.     memcpy (x0_p, x, element_size);
  139.   }
  140.  
  141.   if (copyfunc) {
  142.     destructor(x);
  143.     destructor(new_x);
  144.   } else {
  145.     free (x);
  146.     free (new_x);
  147.   }
  148. }
  149.  
  150. /* implementation of a simulated annealing algorithm with many tries */
  151.  
  152. void 
  153. gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
  154.               gsl_siman_step_t take_step,
  155.               gsl_siman_metric_t distance,
  156.               gsl_siman_print_t print_position,
  157.               size_t element_size,
  158.               gsl_siman_params_t params)
  159. {
  160.   /* the new set of trial points, and their energies and probabilities */
  161.   void *x, *new_x;
  162.   double *energies, *probs, *sum_probs;
  163.   double Ex;            /* energy of the chosen point */
  164.   double T;            /* the temperature */
  165.   int i, done;
  166.   double u;            /* throw the die to choose a new "x" */
  167.   int n_iter;
  168.  
  169.   if (print_position) {
  170.     printf ("#-iter    temperature       position");
  171.     printf ("         delta_pos        energy\n");
  172.   }
  173.  
  174.   x = (void *) malloc (params.n_tries * element_size);
  175.   new_x = (void *) malloc (params.n_tries * element_size);
  176.   energies = (double *) malloc (params.n_tries * sizeof (double));
  177.   probs = (double *) malloc (params.n_tries * sizeof (double));
  178.   sum_probs = (double *) malloc (params.n_tries * sizeof (double));
  179.  
  180.   T = params.t_initial;
  181. /*    memcpy (x, x0_p, element_size); */
  182.   memcpy (x, x0_p, element_size);
  183.   done = 0;
  184.  
  185.   n_iter = 0;
  186.   while (!done)
  187.     {
  188.       Ex = Ef (x);
  189.       for (i = 0; i < params.n_tries - 1; ++i)
  190.     {            /* only go to N_TRIES-2 */
  191.       /* center the new_x[] around x, then pass it to take_step() */
  192.       sum_probs[i] = 0;
  193.       memcpy ((char *)new_x + i * element_size, x, element_size);
  194.       take_step (r, (char *)new_x + i * element_size, params.step_size);
  195.       energies[i] = Ef ((char *)new_x + i * element_size);
  196.       probs[i] = exp (-(energies[i] - Ex) / (params.k * T));
  197.     }
  198.       /* now add in the old value of "x", so it is a contendor */
  199.       memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size);
  200.       energies[params.n_tries - 1] = Ex;
  201.       probs[params.n_tries - 1] = exp (-(energies[i] - Ex) / (params.k * T));
  202.  
  203.       /* now throw biased die to see which new_x[i] we choose */
  204.       sum_probs[0] = probs[0];
  205.       for (i = 1; i < params.n_tries; ++i)
  206.     {
  207.       sum_probs[i] = sum_probs[i - 1] + probs[i];
  208.     }
  209.       u = gsl_rng_uniform (r) * sum_probs[params.n_tries - 1];
  210.       for (i = 0; i < params.n_tries; ++i)
  211.     {
  212.       if (u < sum_probs[i])
  213.         {
  214.           memcpy (x, (char *)new_x + i * element_size, element_size);
  215.           break;
  216.         }
  217.     }
  218.       if (print_position)
  219.     {
  220.       printf ("%5d\t%12g\t", n_iter, T);
  221.       print_position (x);
  222.       printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex);
  223.     }
  224.       T /= params.mu_t;
  225.       ++n_iter;
  226.       if (T < params.t_min)
  227.     {
  228.       done = 1;
  229.     }
  230.     }
  231.  
  232.   /* now return the value via x0_p */
  233.   memcpy (x0_p, x, element_size);
  234.  
  235.   /*  printf("the result is: %g (E=%g)\n", x, Ex); */
  236.  
  237.   free (x);
  238.   free (new_x);
  239.   free (energies);
  240.   free (probs);
  241.   free (sum_probs);
  242. }
  243.